/*
 * =====================================================================================
 *
 *       Filename:  C_T_P3_2D_2D_2D.h
 *
 *    Description:  define the base functions in 2D case
 *
 *        Version:  1.0
 *        Created:  2012/04/03 03时02分41秒
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:   (), @lsec.cc.ac.cn
 *        Company:  
 *
 * =====================================================================================
 */
#ifndef __CTP32D2D__
#define __CTP32D2D__

#include <stdbool.h>
#include "mesh.h"
#include "enumerations.h"
#include "constants.h"
// ***********************************************************************
// (P3)^2 element, conforming, 2D
// ***********************************************************************
static INT C_T_P3_2D_2D_dof[3] = {2,4,2};
static INT C_T_P3_2D_2D_Num_Bas   = 20;
static INT C_T_P3_2D_2D_Value_Dim = 2;
static INT C_T_P3_2D_2D_Inter_Dim = 2;
static INT C_T_P3_2D_2D_Polydeg   = 3;
static bool C_T_P3_2D_2D_IsOnlyDependOnRefCoord = 1;
static INT C_T_P3_2D_2D_Accuracy  = 3;
static MAPTYPE C_T_P3_2D_2D_Maptype = Affine;

static void C_T_P3_2D_2D_InterFun(DOUBLE RefCoord[], INT dim, DOUBLE *values)
{
   DOUBLE x, y;
   x=RefCoord[0]; y=RefCoord[1];
   
     INT i; 
for(i=0;i<40;i++)    values[i] = 0.0;
for(i=0;i<2;i++)
{
  values[4*0+i*3] =18.0*x*y - (11*y)/2 - (11*x)/2 - (27*x*y*y)/2 - (27*x*x*y)/2 + 9*x*x - (9*x*x*x)/2 + 9*y*y - (9*y*y*y)/2 + 1;
   values[4*1+i*3] = x - (9*x*x)/2 + (9*x*x*x)/2;
   values[4*2+i*3] = y - (9*y*y)/2 + (9*y*y*y)/2;
   values[4*3+i*3] = (27*x*x*y)/2 - (9*x*y)/2;
   values[4*4+i*3] = (27*x*y*y)/2 - (9*x*y)/2;
   values[4*5+i*3] = (9*x*y)/2 - (9*y)/2 - (27*x*y*y)/2 + 18*y*y - (27*y*y*y)/2;
   values[4*6+i*3] = 9*y - (45*x*y)/2 + 27*x*y*y + (27*x*x*y)/2 - (45*y*y)/2 + (27*y*y*y)/2;
   values[4*7+i*3] = 9*x - (45*x*y)/2 + (27*x*y*y)/2 + 27*x*x*y - (45*x*x)/2 + (27*x*x*x)/2;
   values[4*8+i*3] = (9*x*y)/2 - (9*x)/2 - (27*x*x*y)/2 + 18*x*x - (27*x*x*x)/2;
   values[4*9+i*3] = 27*x*y - 27*x*y*y - 27*x*x*y;
      }
}


static void C_T_P3_2D_2D_D00(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y;
   x=RefCoord[0]; y=RefCoord[1];
     INT i; 
for(i=0;i<40;i++)    values[i] = 0.0;
for(i=0;i<2;i++)
{
  values[4*0+i*3] =18.0*x*y - (11*y)/2 - (11*x)/2 - (27*x*y*y)/2 - (27*x*x*y)/2 + 9*x*x - (9*x*x*x)/2 + 9*y*y - (9*y*y*y)/2 + 1.0;
   values[4*1+i*3] = x - (9*x*x)/2 + (9*x*x*x)/2;
   values[4*2+i*3] = y - (9*y*y)/2 + (9*y*y*y)/2;
   values[4*3+i*3] = (27*x*x*y)/2 - (9*x*y)/2;
   values[4*4+i*3] = (27*x*y*y)/2 - (9*x*y)/2;
   values[4*5+i*3] = (9*x*y)/2 - (9*y)/2 - (27*x*y*y)/2 + 18*y*y - (27*y*y*y)/2;
   values[4*6+i*3] = 9*y - (45*x*y)/2 + 27*x*y*y + (27*x*x*y)/2 - (45*y*y)/2 + (27*y*y*y)/2;
   values[4*7+i*3] = 9*x - (45*x*y)/2 + (27*x*y*y)/2 + 27*x*x*y - (45*x*x)/2 + (27*x*x*x)/2;
   values[4*8+i*3] = (9*x*y)/2 - (9*x)/2 - (27*x*x*y)/2 + 18*x*x - (27*x*x*x)/2;
   values[4*9+i*3] = 27*x*y - 27*x*y*y - 27*x*x*y;
   }
}
static void C_T_P3_2D_2D_D10(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y;
   x=RefCoord[0]; y=RefCoord[1];
     INT i; 
for(i=0;i<40;i++)    values[i] = 0.0;
for(i=0;i<2;i++)
{
  values[4*0+i*3] =18*x + 18*y - 27*x*y - (27*x*x)/2 - (27*y*y)/2 - 5.5;
   values[4*1+i*3] = (27*x*x)/2 - 9*x + 1.0;
   values[4*2+i*3] = 0.0;
   values[4*3+i*3] = 27*x*y - 4.5*y;
   values[4*4+i*3] = (27*y*y)/2 - 4.5*y;
   values[4*5+i*3] = (9*y)/2 - (27*y*y)/2;
   values[4*6+i*3] = 27*x*y - (45*y)/2 + 27*y*y;
   values[4*7+i*3] = 54*x*y - (45*y)/2 - 45*x   + (81*x*x)/2 + (27*y*y)/2 + 9.0;
   values[4*8+i*3] = 36*x   + (9*y)/2  - 27*x*y - (81*x*x)/2 - 4.5;
   values[4*9+i*3] = 27*y - 54*x*y - 27*y*y;
   }
}
static void C_T_P3_2D_2D_D01(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y;
   x=RefCoord[0]; y=RefCoord[1];
     INT i; 
for(i=0;i<40;i++)    values[i] = 0.0;
for(i=0;i<2;i++)
{
  values[4*0+i*3] =18*x + 18*y - 27*x*y - (27*x*x)/2 - (27*y*y)/2 - 5.5;
   values[4*1+i*3] = 0.0;
   values[4*2+i*3] = (27*y*y)/2 - 9*y + 1.0;
   values[4*3+i*3] = (27*x*x)/2 - (9*x)/2;
   values[4*4+i*3] = 27*x*y - (9*x)/2;
   values[4*5+i*3] = (9*x)/2 + 36*y - 27*x*y - (81*y*y)/2 - 4.5;
   values[4*6+i*3] = 54*x*y - 45*y - (45*x)/2 + (27*x*x)/2 + (81*y*y)/2 + 9.0;
   values[4*7+i*3] = 27*x*y - (45*x)/2 + 27*x*x;
   values[4*8+i*3] = (9*x)/2 - (27*x*x)/2;
   values[4*9+i*3] = 27*x - 54*x*y - 27*x*x;
   }
}
static void C_T_P3_2D_2D_D20(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y;
   x=RefCoord[0]; y=RefCoord[1];
     INT i; 
for(i=0;i<40;i++)    values[i] = 0.0;
for(i=0;i<2;i++)
{
  values[4*0+i*3] =18.0 - 27.0*y - 27.0*x;
   values[4*1+i*3] = 27.0*x - 9.0;
   values[4*2+i*3] = 0.0;
   values[4*3+i*3] = 27*y;
   values[4*4+i*3] = 0.0;
   values[4*5+i*3] = 0.0;
   values[4*6+i*3] = 27*y;
   values[4*7+i*3] = 81*x + 54*y - 45.0;
   values[4*8+i*3] = 36.0 - 27.0*y - 81.0*x;
   values[4*9+i*3] = -54.0*y;
   }
}
static void C_T_P3_2D_2D_D11(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y;
   x=RefCoord[0]; y=RefCoord[1];
     INT i; 
for(i=0;i<40;i++)    values[i] = 0.0;
for(i=0;i<2;i++)
{
  values[4*0+i*3] =18.0 - 27*y - 27*x;
   values[4*1+i*3] = 0.0;
   values[4*2+i*3] = 0.0;
   values[4*3+i*3] = 27*x - 4.5;
   values[4*4+i*3] = 27*y - 4.5;
   values[4*5+i*3] = 4.5 - 27*y;
   values[4*6+i*3] = 27*x + 54*y - 22.5;
   values[4*7+i*3] = 54*x + 27*y - 22.5;
   values[4*8+i*3] = 4.5 - 27*x;
   values[4*9+i*3] = 27.0 - 54*y - 54*x;
   }
}
static void C_T_P3_2D_2D_D02(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
   DOUBLE x, y;
   x=RefCoord[0]; y=RefCoord[1];
     INT i; 
for(i=0;i<40;i++)    values[i] = 0.0;
for(i=0;i<2;i++)
{
  values[4*0+i*3] =18.0 - 27.0*y - 27.0*x;
   values[4*1+i*3] = 0.0;
   values[4*2+i*3] = 27*y - 9.0;
   values[4*3+i*3] = 0.0;
   values[4*4+i*3] = 27*x;
   values[4*5+i*3] = 36.0 - 81.0*y - 27.0*x;
   values[4*6+i*3] = 54*x + 81.0*y - 45.0;
   values[4*7+i*3] = 27.0*x;
   values[4*8+i*3] = 0.0;
   values[4*9+i*3] = -54.0*x;
   }
}
 

static void C_T_P3_2D_2D_Nodal(ELEMENT * elem, FUNCTIONVEC *fun, INT dim, DOUBLE* values)
{ 
    DOUBLE coord[20] =   {0.0, 0.0,  1.0, 0.0,    0.0,1.0, 
                  2.0/3.0,  1.0/3.0,   1.0/3.0,  2.0/3.0,  
                  0.0, 2.0/3.0, 0.0, 1.0/3.0, 
                  1.0/3.0, 0.0,  2.0/3.0, 0.0,
		            1.0/3.0,    1.0/3.0}; 
    INT i;
    DOUBLE w0,w1,w2, COORD[2];
    for(i=0;i<10;i++){
        w1 = coord[2*i];
        w2 = coord[2*i+1];
        w0 = 1.0-w1-w2;
        COORD[0] = w0*elem->Vert_X[0]+w1*elem->Vert_X[1]+w2*elem->Vert_X[2];
        COORD[1] = w0*elem->Vert_Y[0]+w1*elem->Vert_Y[1]+w2*elem->Vert_Y[2];
        fun(COORD,dim,values+i*dim);
    }
}
#endif